#!/usr/bin/perl
#Programmer: Rory Carmichael
#Purpose: Scan min-max results to identify "peaks" (handles for rare codon clusters) and produces a file noting the distance of each window from the nearest peak
#Usage: better-peak-find.pl <inputfile> <cutoff>

use List::Util qw{min max};
if ($#ARGV != 1) {
	die("Usage: better-peak-find.pl <input %min-max file> <cutoff>\n");
}
my $inputfile = shift;
my $cutoff = shift;
if (-e $cutoff) {
	$cutoff = `cat $cutoff`;
}

sub printdists {
	my $curid = shift;
	my $curstring = shift;
	if (defined($curid)) {
		my @peak;
		my @vals = split(/,/, $curstring);
		#the edge condition here kind of takes into account the fact that there will be a trailing , on each line... sorry for the obfuscated math
		#This loop takes care of the peak stuff
		my $leftdist = 0;
#		print STDERR ">INITIAL DATA\n" . join(',',@vals) . "\n";
		for (my $i=0; $i<=$#vals; $i++) {
			my $lower = max(0,$i-5);
			my $upper = min($#vals, $i);
			my @window = @vals[$lower..$upper];
			if (($vals[$i] == min(@window)) && ($vals[$i] < $cutoff)) {
#				print STDERR "SAW A PEAK! $i\n";
#				print STDERR "Value is: $vals[$i]\n";
#				print STDERR "Window is: " . join(',', @window) . "\n";
				push @peak, 0;
				$leftdist=0;
			} else {
				push @peak, ++$leftdist;
			}
		}
#		print STDERR ">INTERMEDIATE DATA\n" . join(',',@peak) . "\n";
		#At this stage the array is populated with distances from the peak to the left (or the left bound of the array if there are no peaks to the left)
		if($#vals+1==$peak[$#vals]) {
			print STDERR "AINT NO PEAKS HERE!\n";
			#No peaks... not sure what to do here
			#EDIT: we set all the peak distances equal to the length of the sequence
			for (my $i=$#vals;$i>=0;$i--) {
				$peak[$i] = $peak[$#vals];
			}
		} else {
			my $numpeaks = 0;
			my $rightdist = 0;
			#This loop takes care of the right distance stuff
			for (my $i=$#vals;$i>=0;$i--) {
				if ($peak[$i] == 0) {
					$numpeaks++;
					$rightdist = 0;
				} elsif (($i+1 == $peak[$i]) || ($numpeaks > 0 && $rightdist < $peak[$i])) {
				#if we are closer to a real peak on the right than on the left, choose the right distance
					$peak[$i] = ++$rightdist;
				}
			}
		}
		print "$curid\n" . join(',', @peak) . "\n";
	}
}

open(IF, $inputfile);
my $curid;
my $curstring;
while (my $line = <IF>) {
	chomp $line;
	if ($line =~ /^>/) {
		print STDERR "looking at $line\n";
		&printdists($curid,$curstring);
		$curid = $line;
		$curstring = "";
	} else {
		$curstring .= $line;
	}
}
&printdists($curid,$curstring);
close(IF)
